Part 2 - ASA Short Course
Second-generation p-values (Full-Day)
Introduction
- Jeffrey D. Blume, PhD
- School of Data Science, University of Virginia
- Megan H. Murray, PhD
- Department of Biostatistics, Vanderbilt University
Resources:
- GitHub with Slides and Code: www.github.com/murraymegan/SGPV-ASA-Short-Course
- RStudio Desktop: www.rstudio.com/products/rstudio/download
- Interrupt or use Zoom chat for questions!
- For technical difficulties email Megan: megan.c.hollister@vanderbilt.edu
SGPV Package
Functions Included:
sgpvalue()- This function computes the second-generation p-value (SGPV) and its associated delta gaps, as introduced in Blume et al. (2018).
sgpower()- Calculate power and type I error values from significance testing based on second-generation p-values as the inferential metric.
plotsgpv()- This function displays user supplied interval estimates (support intervals, confidence intervals, credible intervals, etc.) according to its associated second-generation p-value ranking.
plotman()- This function displays a modified Manhattan-style plot colored according to second-generation p-value status.
fdrisk()- This function computes the false discovery risk (sometimes called the “empirical bayes FDR”) for a second-generation p-value of 0, or the false confirmation risk for a second-generation p-value of 1.
1 Part 2
1.1 Part 2a
1.1.1 Statistical Properties
3 inference outcomes:
- Data are consistent with null \(p_\delta=1\)
- Data are consistent with alternative \(p_\delta=0\)
- Data are inconclusive \(0<p_\delta<1\)
2 underlying truths:
- Null is true
- Alternative is true
prob.alt.sg = function(n, null.delta,
theta0=0,
theta=0,
V=1,
alpha=0.05){
results = rep(NA, length(n))
results = pnorm(sqrt(n)*(theta0-null.delta)/sqrt(V)-sqrt(n)*theta/sqrt(V)-qnorm(1-alpha/2))+
pnorm(-1*sqrt(n)*(theta0+null.delta)/sqrt(V)+sqrt(n)*theta/sqrt(V)-qnorm(1-alpha/2))
results
}prob.null.sg = function(n,
null.delta,
theta0=0,
theta=0,
V=1,
alpha=0.05){
results = rep(NA, length(n))
if(length(null.delta)>1){
con.large = which(null.delta>(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
con.small = which(null.delta<=(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
results[con.small] = rep(0, length(con.small))
results[con.large] =
pnorm(sqrt(n[con.large])*theta0+sqrt(n[con.large])*null.delta[con.large]/sqrt(V)-sqrt(n[con.large])*theta/sqrt(V)-qnorm(1-alpha/2))-
pnorm(sqrt(n[con.large])*theta0-sqrt(n[con.large])*null.delta[con.large]/sqrt(V)-sqrt(n[con.large])*theta/sqrt(V)+qnorm(1-alpha/2))
}else if(null.delta>(qnorm(1-alpha/2)*sqrt(V)/sqrt(n))){
results =
pnorm(sqrt(n)*theta0+sqrt(n)*null.delta/sqrt(V)-sqrt(n)*theta/sqrt(V)-qnorm(1-alpha/2))-
pnorm(sqrt(n)*theta0-sqrt(n)*null.delta/sqrt(V)-sqrt(n)*theta/sqrt(V)+qnorm(1-alpha/2))
}else{
results = rep(0, length(theta))
}
results
}prob.incon.sg = function(n, null.delta,
theta0=0,
theta=0,
V=1,
alpha=0.05){
results = rep(NA, length(n))
con.large = which(null.delta>(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
con.small = which(null.delta<=(qnorm(1-alpha/2)*sqrt(V)/sqrt(n)))
results = 1-(prob.alt.sg(n, null.delta,
theta0,
theta,
V,
alpha)+
prob.null.sg(n, null.delta,
theta0,
theta,
V,
alpha))
results
}theta.1 = seq(-3,3,by=0.01)
plot(theta.1, prob.alt.sg(n=50,
null.delta = 0,
theta=theta.1), type="l",
ylim=c(0,1),
xlab="Alternative",
ylab="Probability",
main="Prob of being consistent with alternative")
lines(theta.1, prob.alt.sg(n=50,
null.delta = 0.5,
theta=theta.1),
col="green")
lines(theta.1, prob.alt.sg(n=50,
null.delta = 1,
theta=theta.1),
col="blue")
abline(h=0.05, lty=2)
legend("bottomright",
legend=c("Delta=0", "Delta=0.5", "Delta=1"),
col=c("black", "green","blue"), lty=1)
theta.1 = seq(-3,3,by=0.01)
plot(theta.1, prob.null.sg(n=50,
null.delta = 0,
theta=theta.1), type="l",
ylim=c(0,1),
xlab="Alternative",
ylab="Probability",
main="Prob of being consistent with null")
lines(theta.1, prob.null.sg(n=50,
null.delta = 0.5,
theta=theta.1),
col="green")
lines(theta.1, prob.null.sg(n=50,
null.delta = 1,
theta=theta.1),
col="blue")
abline(h=0.05, lty=2)
legend("bottomright",
legend=c("Delta=0", "Delta=0.5", "Delta=1"),
col=c("black", "green","blue"), lty=1)
theta.1 = seq(-3,3,by=0.01)
plot(theta.1, prob.incon.sg(n=50,
null.delta = 0,
theta=theta.1), type="l",
ylim=c(0,1),
xlab="Alternative",
ylab="Probability",
main="Prob of being inconclusive")
lines(theta.1, prob.incon.sg(n=50,
null.delta = 0.5,
theta=theta.1),
col="green")
lines(theta.1, prob.incon.sg(n=50,
null.delta = 1,
theta=theta.1),
col="blue")
abline(h=0.05, lty=2)
legend("bottomright",
legend=c("Delta=0", "Delta=0.5", "Delta=1"),
col=c("black", "green","blue"), lty=1)
1.2 Part 2b
1.2.1 Equivalence tests: TOST
- Establish bioequivalence between data and an established range
- Example: A pharmaceutical company tests for drug approval by comparing new drug’s performance to an approved drug’s performance
- Tests are ordinary, one-sided, α-level t-tests
- Flips the null and alternative hypotheses
n=6
dat = rnorm(n, 0, 1)
a = t.test(dat)$conf.int[1]
b = t.test(dat)$conf.int[2]
x_bar = mean(dat)
S = sd(dat)
Z = (b-x_bar)/(S/sqrt(n))
theta_p = 0.5*0.75
theta_m = -0.5*0.75
out_TOST = tsum_TOST(m1=x_bar,
mu=0,
sd1=S,
n1=n,
low_eqbound=theta_m,
high_eqbound=theta_p,
alpha=0.05,
eqbound_type = "raw")
out_TOST$TOST## t SE df p.value
## t-test -0.4206392 0.1519841 5 0.69148580
## TOST Lower 2.0467240 0.1519841 5 0.04802133
## TOST Upper -2.8880023 0.1519841 5 0.01713328max(out_TOST$TOST$p.value[2], out_TOST$TOST$p.value[3])## [1] 0.04802133out_sgpv = sgpvalue(est.lo=a,
est.hi=b,
null.lo=theta_m,
null.hi=theta_p)
out_sgpv ## $p.delta
## [1] 0.8981052
##
## $delta.gap
## [1] NAset.seed(999)
iter=500
n=6
results = as.data.frame(matrix(ncol=3, nrow=iter))
colnames(results) = c("sgpv", "tost", "Case")
half.case = NULL
half.case2 = NULL
for(i in 1:iter){
dat = rnorm(n, 0, 1)
a = t.test(dat)$conf.int[1]
b = t.test(dat)$conf.int[2]
x_bar = mean(dat)
S = sd(dat)
Z = (b-x_bar)/(S/sqrt(n))
theta_p = 0.5*0.75
theta_m = -0.5*0.75
out_TOST = tsum_TOST(m1=x_bar,
mu=0,
sd1=S,
n1=n,
low_eqbound=theta_m,
high_eqbound=theta_p,
alpha=0.05,
eqbound_type = "raw")
out_sgpv = sgpvalue(est.lo=a,
est.hi=b,
null.lo=theta_m,
null.hi=theta_p)$p.delta
results[i,1] = out_sgpv
results[i,2] = pmax(out_TOST$TOST$p.value[2],out_TOST$TOST$p.value[3])
if(theta_m<a & theta_p>a & theta_p<b){
# Case 4 some overlap
results[i,3] = "\n Case 4: \nPartial Overlap"
}else if(theta_m<b & theta_p>b & theta_m>a){
# Case 3 some overlap
results[i,3] = "\n Case 3: \nPartial Overlap"
}else if(theta_m>b | theta_p<a){
# Case 1 no overlap
results[i,3] = "\n Case 1: \nNo Overlap"
}else if((theta_m>a & theta_p<b)|(theta_m<=a & theta_p>=b)){
# Case 2 complete overlap
results[i,3] = "\n Case 2: \nComplete Overlap"
}
}
library(ggplot2)
scaleFUN <- function(x) sprintf("%.1f", x)
g = ggplot(results, aes(sgpv, tost))+
xlab("SGPVs") +
ylab("Equivalence p-values") +
geom_point(aes(alpha=0.2)) +
geom_abline(intercept = 1, slope = -1,linetype='dashed', alpha=0.6)+
xlim(c(0,1))+
scale_x_continuous(breaks= c( 0,0.120001,0.880001,1),labels=c("0", "~0.1", "~0.9", "1"), minor_breaks = c(0.12,0.88),
sec.axis=sec_axis(~./1, name="",breaks= c( 0.04,0.51,0.98), labels=c("Consistent \n with Alternative", "\n Inconclusive", "Consistent \n with Null")))+
scale_y_continuous(breaks= c(0,0.07,1), labels=c("0", expression(alpha), "1"),minor_breaks = c(0.070001),
sec.axis=sec_axis(~./1, name="", breaks= c(0.005,0.61), labels=c("Consistent \n with Null","Inconclusive")))+
annotate(geom="text", x=0, y=0.9, label="D",
size=5,
family = "Courier-Bold")+
annotate(geom="text", x=0.85, y=0.98, label="E",
size=5,
family = "Courier-Bold")+
annotate(geom="text", x=0.91, y=0.98, label="F",
size=5,
family = "Courier-Bold")+
annotate(geom="text", x=0, y=0, label="H",
size=5,
family = "Courier-Bold")+
annotate(geom="text", x=0.15, y=0, label="I",
size=5,
family = "Courier-Bold")+
annotate(geom="text", x=0.91, y=0, label="J",
size=5,
family = "Courier-Bold")+
theme_bw()+
theme(axis.text.y = element_text(size=10, hjust=-0.5),
axis.text.y.right = element_text(size=10, hjust=0.7),
axis.text.x.top = element_text(size=10, vjust=2),
axis.text.x = element_text(size=10, vjust=-0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_line(size = 0.8,colour="darkgrey"))+
guides(alpha=FALSE)
g
ggplot(results, aes(sgpv, tost, col=Case, pch=Case))+
geom_point(aes(), cex=1.5) +
# geom_jitter(width = 0.02, height = 0.02)+
# ggtitle(paste0("Derived Sample Size= ",n," Iterations= ", iter))+
xlab("SGPV") +
ylab("TOST Eq PV") +
# facet_wrap(~ Case) +
geom_abline(intercept = 1, slope = -1,linetype='dashed', alpha=0.6)+
scale_color_manual(values=c("#0000B9",
"#00A300",
"#FD49A0",
"#FFA500",
"black"))+
scale_shape_manual(values=c(16,21,17,24,21))+
guides(alpha=FALSE)+
scale_y_continuous(labels=scaleFUN)+
scale_x_continuous(labels=scaleFUN)+
guides(color = guide_legend(byrow = TRUE))+
theme_bw()+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.spacing.y = unit(0.3, 'cm'))
1.3 Part 2c
1.3.1 SGPV False Discovery Risk
- False Discovery Rate (FDR)
- SGPV = 0
- False Confirmation Rate (FCR)
- SGPV = 1
sgpv::fdrisk()- This function computes the false discovery risk (sometimes called the “empirical bayes FDR”) for a second-generation p-value of 0, or the false confirmation risk for a second-generation p-value of 1.
###### FDR rates
# false discovery risk with 95% confidence level
fdrisk(sgpval = 0,
null.lo = log(1/1.1),
null.hi = log(1.1),
std.err = 0.8,
null.weights = 'Uniform',
null.space = c(log(1/1.1), log(1.1)),
alt.weights = 'Uniform',
alt.space = 2 + c(-1,1)*qnorm(1-0.05/2)*0.8,
interval.type = 'confidence',
interval.level = 0.05)## [1] 0.05949861## with truncated normal weighting distribution
fdrisk(sgpval = 0,
null.lo = log(1/1.1),
null.hi = log(1.1),
std.err = 0.8,
null.weights = 'Point',
null.space = 0,
alt.weights = 'TruncNormal',
alt.space = 2 + c(-1,1)*qnorm(1-0.041/2)*0.8,
interval.type = 'likelihood', interval.level = 1/8)## [1] 0.04902575# false discovery risk with LSI and wider null hypothesis
fdrisk(sgpval = 0,
null.lo = log(1/1.5),
null.hi = log(1.5),
std.err = 0.8,
null.weights = 'Point',
null.space = 0,
alt.weights = 'Uniform',
alt.space = 2.5 + c(-1,1)*qnorm(1-0.041/2)*0.8,
interval.type = 'likelihood', interval.level = 1/8)## [1] 0.01688343# false confirmation risk example
fdrisk(sgpval = 1,
null.lo = log(1/1.5),
null.hi = log(1.5),
std.err = 0.15,
null.weights = 'Uniform',
null.space = 0.01 + c(-1,1)*qnorm(1-0.041/2)*0.15,
alt.weights = 'Uniform',
alt.space = c(log(1.5), 1.25*log(1.5)),
interval.type = 'likelihood', interval.level = 1/8)## [1] 0.03059522###
##
#False Discovery and Confirmation Rates
power.pv=function(n=20,sd=1,alpha=0.05,mu.l=-3,mu.u=3,delta=0.3,mu.0=0,r=1){
## Computes power when delta is zero and non-zero
## Inputs:
## n is sample size
## indifference zone is: mu.0 +/- delta
## alternative space is: mu in [mu.l, mu.u]
## data ~ N(mu,sd)
## r = P(H1)/P(H0) for simple H1 and H0
mu=seq(mu.l,mu.u,0.01)
## Power (null is mu.0 ; alt is mu)
pow1.d=pnorm(-sqrt(n)*(delta-(mu.0-mu))/sd-qnorm(1-alpha/2))
pow2.d=pnorm(-sqrt(n)*(delta+(mu.0-mu))/sd-qnorm(1-alpha/2))
power.d=pow1.d+pow2.d
## Power without indifference zone (delta=0)
pow1.0=pnorm(-sqrt(n)*(0-(mu.0-mu))/sd-qnorm(1-alpha/2))
pow2.0=pnorm(-sqrt(n)*(0+(mu.0-mu))/sd-qnorm(1-alpha/2))
power.0=pow1.0+pow2.0
## Type I error (null is mu.0 ; alt is mu.0)
alpha.d=2*pnorm(-sqrt(n)*(delta)/sd-qnorm(1-alpha/2))
## Type I error without indifference zone (delta=0)
alpha.0=2*pnorm(-sqrt(n)*(0)/sd-qnorm(1-alpha/2))
## P(sgpv=1|H0) not quite alpha b/c rules out inconclusive
gam1.0=pnorm(sqrt(n)*(delta)/sd-qnorm(1-alpha/2))
gam2.0=pnorm(sqrt(n)*(-delta)/sd+qnorm(1-alpha/2))
gam.0=(gam1.0-gam2.0)*(delta>qnorm(1-alpha/2)*sd/sqrt(n))
##
gam1.a=pnorm(sqrt(n)*(delta+mu.0-mu)/sd-qnorm(1-alpha/2))
gam2.a=pnorm(sqrt(n)*(mu.0-delta-mu)/sd+qnorm(1-alpha/2))
gam.a=(gam1.a-gam2.a)*(delta>qnorm(1-alpha/2)*sd/sqrt(n))
## False discovery rates (fdr1=FDR ; fdr0=FNDR)
fdr.d = 1 / (1 + power.d/alpha.d * r)
fndr.0 = 1 / (1 + (1-alpha.0)/(1-power.0) * (1/r))
fdr.0 = 1 / (1 + power.0/alpha.0 * r)
fndr.d = 1 / (1 + gam.0/gam.a * (1/r))
ret=cbind(mu,alpha.d,power.d,fdr.d,fndr.d,
alpha.0,power.0,fdr.0,fndr.0,gam.0,gam.a)
colnames(ret)=c('mu','alpha.d','power.d','fdr.d','fndr.d',
'alpha.0','power.0','fdr.0','fndr.0',
'gam.0','gam.a')
ret
}par(mfrow=c(1,1))
prev.ratio=1
set.delta=0.5
set.n=16
all=power.pv(delta=set.delta,n=set.n,r=prev.ratio,alpha=0.05)
plot(all[,'mu'],all[,'fdr.0'],type="n",ylab="Probability",xlab="Alternative in SD units",
ylim=c(0,prev.ratio/(prev.ratio+1)),xlim=c(-3,3))
title(main="False discovery and confirmation rates")
lines(all[,'mu'],all[,'fdr.0'],col="firebrick",lwd=1,lty=2)
lines(all[,'mu'],all[,'fndr.0'],col="dodgerblue",lwd=1,lty=2)
lines(all[,'mu'],all[,'fdr.d'],col="firebrick",lwd=2)
lines(all[,'mu'],all[,'fndr.d'],col="dodgerblue",lwd=2)
legend("topleft",c("",
expression(paste("P-value < ",0.05,sep="")),
expression(paste(P[delta]," = 0",sep="")),
"", "",
expression(paste("P-value > ",0.05,sep="")),
expression(paste(P[delta]," = 1",sep=""))),
col=c("","firebrick","firebrick","","",
"dodgerblue","dodgerblue"),
lty=c(NA,2,1,NA,NA,2,1),lwd=c(NA,1,2,NA,NA,1,2),bty="n")
text(-2.875,0.5,"FDR")
text(-2.875,0.42,"FCR")
axis(side = 1, at = c(-set.delta,set.delta),cex.axis=0.8,col.ticks="black",lwd.ticks=1.5,tcl=0.6,labels=FALSE)
axis(side = 3, at = c(-set.delta,set.delta),cex.axis=0.8,col.ticks="black",lwd.ticks=1.5,tcl=0.6,labels=FALSE)
text(0,-0.0125,expression(paste(H[0]," Zone",sep="")),cex=0.8)
text(0,prev.ratio/(prev.ratio+1)*1.015,expression(paste(H[0]," Zone",sep="")),cex=0.8)
text(2.5,0.05/(prev.ratio+0.05)+.01,"alpha/(alpha+r)")
1.3.2 FDRestimation
Download the package from:
Our corresponding paper: https://f1000research.com/articles/10-441
p.fdr.output = p.fdr(leukstats$p.value)
head(p.fdr.output$`Results Matrix`)## BH FDRs Adjusted p-values Raw p-values
## 1 0.1892484 0.1892484 0.05615324
## 2 0.9762477 0.9759186 0.93296848
## 3 0.6970430 0.6969894 0.47936375
## 4 0.2141365 0.2141365 0.06684255
## 5 0.9160462 0.9156422 0.81439179
## 6 0.6859496 0.6859383 0.46673059summary(p.fdr.output)##
## Call:
## p.fdr(pvalues = leukstats$p.value)
##
## Number of tests: 7128
## Raw p-value Range: [0, 0.99964]
##
## Adjustment Method: BH
## False Discovery Rate Range: [0, 1]
## Findings at 0.05 level:
## Significant (Reject): 1233
## Inconclusive (Fail to Reject): 5895
## ---
## Estimated/Assumed Null Proportion (pi0): 1Plots
plot(p.fdr.output)
plot(p.fdr.output, xlim=c(1100,1600), ylim=c(0, 0.2))
which(-1*(p.fdr.output$`Results Matrix`$`Adjusted p-values`)+p.fdr.output$fdrs>0.001)## integer(0)#Benjamini-Yeukatelli
p.fdr.output = p.fdr(leukstats$p.value, adjust.method = "BY")
head(p.fdr.output$`Results Matrix`)## BY FDRs BY Adjusted p-values Raw p-values BH FDRs
## 1 0.1892484 0.1892484 0.05615324 0.1892484
## 2 1.0000000 1.0000000 0.93296848 0.9762477
## 3 1.0000000 1.0000000 0.47936375 0.6970430
## 4 0.4461177 0.4461177 0.06684255 0.2141365
## 5 1.0000000 1.0000000 0.81439179 0.9160462
## 6 1.0000000 1.0000000 0.46673059 0.6859496plot(p.fdr.output)
plot(p.fdr.output, xlim=c(2000,2400), legend.on = FALSE)
Simple Example
pvalues=c(0.005,0.049,0.05,0.051,0.7)
zvalues=qnorm(pvalues/2, lower.tail = FALSE)
p.fdr.output = p.fdr(pvalues)
adj.pvalues= p.fdr.output$`Results Matrix`$`Adjusted p-values`
adj.fdrs= p.fdr.output$fdrs
single.fdr = c(p.fdr(pvalues=pvalues[1],zvalue=zvalues[1]),
p.fdr(pvalues=pvalues[2],zvalue=zvalues[2]),
p.fdr(pvalues=pvalues[3],zvalue=zvalues[3]),
p.fdr(pvalues=pvalues[4],zvalue=zvalues[4]),
p.fdr(pvalues=pvalues[5],zvalue=zvalues[5]))df = data.frame("Raw p-values"= pvalues,
"Z-values" = zvalues,
"BH Adj p-values" = adj.pvalues,
"BH FDRs" = adj.fdrs,
"Single lower bound FDRs" = single.fdr)
colnames(df) = c("Raw p-values",
"Z-values",
"BH Adj p-values" ,
"BH FDRs" ,
"Single lower bound FDRs" )
kable(round(df,3))%>%
kable_styling(c("striped", "bordered"))| Raw p-values | Z-values | BH Adj p-values | BH FDRs | Single lower bound FDRs |
|---|---|---|---|---|
| 0.005 | 2.807 | 0.025 | 0.025 | 0.019 |
| 0.049 | 1.969 | 0.064 | 0.122 | 0.126 |
| 0.050 | 1.960 | 0.064 | 0.083 | 0.128 |
| 0.051 | 1.951 | 0.064 | 0.064 | 0.130 |
| 0.700 | 0.385 | 0.700 | 0.700 | 0.481 |
Compare to p.adjust
set.seed(999)
pi0 <- 0.8
pi1 <- 1-pi0
n <- 10
n.0 <- ceiling(n*pi0)
n.1 <- n-n.0
sim.data <- c(rnorm(n.1,2,1),rnorm(n.0,0,1))
sim.data.p <- 2*pnorm(-abs(sim.data))
p.adjust.output = p.adjust(sim.data.p, method="fdr")
fdr.output = p.fdr(pvalues=sim.data.p, adjust.method="BH")
head(data.frame("Raw p-values"= sim.data.p,"p.adjust FDRs"=p.adjust.output,"p.fdr adj p-values"=fdr.output$`Results Matrix`$`Adjusted p-values`,"p.fdr FDRs"=fdr.output$fdrs))## Raw.p.values p.adjust.FDRs p.fdr.adj.p.values p.fdr.FDRs
## 1 0.08574923 0.4287462 0.4287462 0.4287462
## 2 0.49180527 0.7025790 0.7025790 0.7025790
## 3 0.42650649 0.7025790 0.7025790 0.7108441
## 4 0.78710602 0.7871060 0.7871060 0.7871060
## 5 0.78154483 0.7871060 0.7871060 0.8683831
## 6 0.57137764 0.7142221 0.7142221 0.7142221plot(rank(sim.data.p, ties="random"),
fdr.output$`Results Matrix`$`Adjusted p-values`,
pch=17,
col="dodgerblue",
cex=2,
ylim=c(0,1),
main="Comparison Plot",
xlab="Rank of p-values",
ylab="")
points(rank(sim.data.p, ties="random"),p.adjust.output, pch=20, col="pink")
points(rank(sim.data.p, ties="random"),fdr.output$fdrs, pch=20, col="firebrick")
legend("bottomright", c("p.adjust FDRs", "p.fdr FDRs","p.fdr Adjusted p-values"), col=c("pink", "firebrick","dodgerblue"), pch=c(20,20,17))
Different Null proportions
get.pi0(leukstats$p.value, estim.method = "last.hist")## [1] 0.6313131get.pi0(leukstats$p.value, estim.method = "storey")## [1] 0.6020751get.pi0(leukstats$p.value, estim.method = "set.pi0", set.pi0=0.8)## [1] 0.8set.seed(999)
pi0 <- 0.8
pi1 <- 1-pi0
n <- 100
n.0 <- ceiling(n*pi0)
n.1 <- n-n.0
sim.data <- c(rnorm(n.1,2,1),rnorm(n.0,0,1))
sim.data.p <- 2*pnorm(-abs(sim.data))
fdr.output = p.fdr(pvalues=sim.data.p, adjust.method="BH")
plot(fdr.output)
fdr.output = p.fdr(pvalues=sim.data.p, adjust.method="BH", set.pi0=0.8)
plot(fdr.output)
get.pi0(sim.data.p, estim.method = "set.pi0", set.pi0=0.7)## [1] 0.7get.pi0(sim.data.p, estim.method = "last.hist")## [1] 1get.pi0(sim.data.p, estim.method = "storey")## [1] 1pi0 <- 0.8
pi1 <- 1-pi0
n <- 10000
n.0 <- ceiling(n*pi0)
n.1 <- n-n.0
sim.data <- c(rnorm(n.1,2,1),rnorm(n.0,0,1))
sim.data.p <- 2*pnorm(-abs(sim.data))
get.pi0(sim.data.p, estim.method = "set.pi0", set.pi0=0.7)## [1] 0.7get.pi0(sim.data.p, estim.method = "last.hist")## [1] 0.796get.pi0(sim.data.p, estim.method = "storey")## [1] 0.8051147Internal Code: Last Histogram Height Method
try.hist = hist(sim.data.p, breaks="scott")
try.mids = try.hist$mids
try.count = try.hist$counts
if(tail(try.mids,1)<0.5){
pi0.hat=0
}else{
pi0.hat = min(tail(try.count,1)*length(try.mids)/sum(try.count),1)
}
pi0.hat## [1] 0.7961.4 Part 2d
1.4.1 Study Planning
Purpose: Evaluate how different techniques of setting the indifference zone effect the SGPV study properties.
What options does a collaborator have when they are uncertain of the indifference zone?
Fixed Intervals
n= seq(1,250, by=0.1)
V=1
alpha=0.05
null.delta0 = rep(0.5, length(n))
null.delta1 = rep(0, length(n))
null.delta2 = rep(1.25, length(n))
null.delta3 = rep(1, length(n))
null.delta4 = rep(0.75, length(n))
null.delta5 = rep(0.25, length(n))
ci.bound = qnorm(1-alpha/2)*sqrt(V)/sqrt(n)plot(n, null.delta1,
type="l",
col="red",
lwd=2,
ylim=c(0,3),
ylab=TeX(r'(Half Interval Width $ (\delta)$)'),
xlab = "Sample Size (n)",
main="")
lines(n, null.delta2,
col="blue",
lwd=2)
lines(n, null.delta3,
col="purple",
lwd=2)
lines(n, null.delta4,
col="darkgreen",
lwd=2)
lines(n, null.delta5,
col="orange",
lwd=2)
lines(n, null.delta0,
col="hotpink",
lwd=2)
lines(n, ci.bound,
col="black",
lty=2,
lwd=2)
# abline(h=0.05, col="black", lty=2)
legend(130,
3,
col=c("black",
"red",
"blue",
"purple",
"darkgreen",
"hotpink",
"orange"),
box.lty=0,
lty=c(2,
1,
1,
1,
1,
1,
1),
cex=0.65,
lwd=2,
legend=c(TeX(r'(Uncertainty Interval $\left(\sqrt{1/n}\right)$)'),
"Point Null: Constant at 0",
"SGPV: Constant at 1.25",
"SGPV: Constant at 1",
"SGPV: Constant at 0.75",
"SGPV: Constant at 0.5",
"SGPV: Constant at 0.25"))
par(mfrow=c(1,2))
par(mar=c(2,2,2,2))
plot(n, prob.alt.sg(n, null.delta1,
theta=0),
type="l",
col="red",
lwd=2,
ylim=c(0,0.055),
ylab="",
main=TeX(r'(Type 1 Error $\beta_0$)'))
lines(n, prob.alt.sg(n, null.delta2,
theta=0),
col="blue",
lwd=2)
lines(n, prob.alt.sg(n, null.delta3,
theta=0),
col="purple",
lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
theta=0),
col="darkgreen",
lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
theta=0),
col="orange",
lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
theta=0),
col="hotpink",
lwd=2)
abline(h=0.05, col="black", lty=2)
plot(n, prob.alt.sg(n, null.delta1,
theta=1),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Power $\beta_1$)'))
lines(n, prob.alt.sg(n, null.delta2,
theta=1),
col="blue",
lwd=2)
lines(n, prob.alt.sg(n, null.delta3,
theta=1),
col="purple",
lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
theta=1),
col="darkgreen",
lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
theta=1),
col="orange",
lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
theta=1),
col="hotpink",
lwd=2)
legend(70, 0.5,
col=c("red",
"blue",
"purple",
"darkgreen",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1,
1),
cex=0.7,
lwd=2,
legend=c("Point Null: Constant at 0",
"SGPV: Constant at 1.25",
"SGPV: Constant at 1",
"SGPV: Constant at 0.75",
"SGPV: Constant at 0.5",
"SGPV: Constant at 0.25"))
par(mfrow=c(1,1))
par(mar=c(4,4,2,2))
plot(n, prob.null.sg(n, null.delta1,
theta=0)+prob.alt.sg(n, null.delta1,
theta=1),
type="l",
col="red",
lwd=2,
ylim=c(0,2),
ylab=TeX(r'( $\beta_1$ + $\omega_0$)'),
main= TeX(r'(Sum of Power and Probability of Null ($\beta_1$ +$\omega_0$))'))
lines(n, prob.null.sg(n, null.delta2,
theta=0)+prob.alt.sg(n, null.delta2,
theta=1),
col="blue",
lwd=2)
# lines(n, prob.null.sg(n, null.delta3,
# theta=0),
# col="purple",
# lwd=2)
lines(n, prob.null.sg(n, null.delta4,
theta=0)+prob.alt.sg(n, null.delta4,
theta=1),
col="darkgreen",
lwd=2)
lines(n, prob.null.sg(n, null.delta5,
theta=0)+prob.alt.sg(n, null.delta5,
theta=1),
col="orange",
lwd=2)
lines(n, prob.null.sg(n, null.delta0,
theta=0)+prob.alt.sg(n, null.delta0,
theta=1),
col="hotpink",
lwd=2)
legend(172, 0.7,
col=c("red",
"blue",
"purple",
"darkgreen",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2),
cex=0.7,
lwd=2,
legend=c(TeX(r'( Constant at 0)'),
TeX(r'(Constant at 1.25)'),
TeX(r'( Constant at 1)'),
TeX(r'( Constant at 0.75)'),
TeX(r'(Constant at 0.5)'),
TeX(r'( Constant at 0.25)')))
par(mfrow=c(1,2))
par(mar=c(2,2,4,2))
plot(n, prob.null.sg(n, null.delta1,
theta=0),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main= TeX(r'(Probability of Null $\omega_0$)'))
lines(n, prob.null.sg(n, null.delta2,
theta=0),
col="blue",
lwd=2)
lines(n, prob.null.sg(n, null.delta3,
theta=0),
col="purple",
lwd=2)
lines(n, prob.null.sg(n, null.delta4,
theta=0),
col="darkgreen",
lwd=2)
lines(n, prob.null.sg(n, null.delta5,
theta=0),
col="orange",
lwd=2)
lines(n, prob.null.sg(n, null.delta0,
theta=0),
col="hotpink",
lwd=2)
plot(n, prob.null.sg(n, null.delta1,
theta=1),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Probability of Null $\omega_1$)'))
lines(n, prob.null.sg(n, null.delta2,
theta=1),
col="blue",
lwd=2)
lines(n, prob.null.sg(n, null.delta3,
theta=1),
col="purple",
lwd=2)
lines(n, prob.null.sg(n, null.delta4,
theta=1),
col="darkgreen",
lwd=2)
lines(n, prob.null.sg(n, null.delta5,
theta=1),
col="orange",
lwd=2)
lines(n, prob.null.sg(n, null.delta0,
theta=1),
col="hotpink",
lwd=2)
legend(70, 0.5,
col=c("red",
"blue",
"purple",
"darkgreen",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1,
1),
cex=0.7,
lwd=2,
legend=c("Point Null: Constant at 0",
"SGPV: Constant at 1.25",
"SGPV: Constant at 1",
"SGPV: Constant at 0.75",
"SGPV: Constant at 0.5",
"SGPV: Constant at 0.25"))
par(mfrow=c(1,2))
par(mar=c(2,2,4,2))
plot(n,prob.incon.sg(n, null.delta1,
theta=0),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Probability of Inconclusive $\gamma_0$)'))
lines(n,prob.incon.sg(n, null.delta2,
theta=0),
col="blue",
lwd=2)
lines(n,prob.incon.sg(n, null.delta3,
theta=0),
col="purple",
lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
theta=0),
col="darkgreen",
lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
theta=0),
col="orange",
lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
theta=0),
col="hotpink",
lwd=2)
plot(n,prob.incon.sg(n, null.delta1,
theta=1),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Probability of Inconclusive $\gamma_1$)'))
lines(n,prob.incon.sg(n, null.delta2,
theta=1),
col="blue",
lwd=2)
lines(n,prob.incon.sg(n, null.delta3,
theta=1),
col="purple",
lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
theta=1),
col="darkgreen",
lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
theta=1),
col="orange",
lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
theta=1),
col="hotpink",
lwd=2)
legend(90, 0.8,
col=c("red",
"blue",
"purple",
"darkgreen",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1,
1),
cex=0.7,
lwd=2,
legend=c("Point Null: Constant at 0",
"SGPV: Constant at 1.25",
"SGPV: Constant at 1",
"SGPV: Constant at 0.75",
"SGPV: Constant at 0.5",
"SGPV: Constant at 0.25"))
Shrinking from 0.9 to 0.45
Let’s see how these look on a plot:
n2= seq(-10,300, by=0.1)
plot(n, ci.bound,
type="n",
col="black",
lty=2,
lwd=2,
ylim=c(-2,2.2),
xlim=c(0,225),
ylab="",
xlab = "Sample Size (n)",
main="",
yaxt = "n")
# lines(n, -1*ci.bound,
# col="black",
# lty=2,
# lwd=2)
lines(n, null.delta2,
col="blue",
lwd=2)
lines(n, -1*null.delta2,
col="blue",
lwd=2)
lines(n, null.delta4,
col="hotpink",
lwd=2)
lines(n, -1*null.delta4,
col="hotpink",
lwd=2)
lines(n, null.delta5,
col="orange",
lwd=2)
lines(n, -1*null.delta5,
col="orange",
lwd=2)
abline(h=k, col="black", lty=1,lwd=2)
abline(h=-1*k, col="black", lty=1,lwd=2)
abline(h=0, col="red", lty=1,lwd=2)
abline(h=1, col="red", lty=2,lwd=2)
polygon(c(n2, rev(n2)), c(rep(-1*k, length(n2)), rev(rep(-1*(k/2), length(n2)))),
border = NA,
col = rgb(210/633,210/633,210/633, alpha = 0.2))
polygon(c(n2, rev(n2)), c(rep((k/2), length(n2)), rev(rep(k, length(n2)))),
border = NA,
col = rgb(210/633,210/633,210/633, alpha = 0.2))
axis(2, at = c(-2,-1,-0.9,-0.45,0,
2,1,0.9,0.45),
labels = c(-2,-1,-0.9,-0.45,TeX(r'($\theta_0=0$)'),
2,TeX(r'($\theta_1=1$)'),0.9,0.45),
las=1)
legend(135,
2.25,
col=c("red",
"red",
"black",
"blue",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
2,
1,
1,
1,
1),
cex=0.9,
lwd=2,
legend=c(TeX(r'(Null Hypothesis $\theta_0=0$)'),
TeX(r'(Alternative Hypothesis $\theta_1=1$)'),
"Original SGPV: Constant at 0.9",
TeX(r'($delta_1(n)\propto n^{-1/4}$)'),
TeX(r'($delta_2(n)\propto n^{-1/3}$)'),
TeX(r'($delta_3(n)\propto n^{-1/2}$)')))
1.4.2 Probability Plots
par(mfrow=c(1,2))
par(mar=c(2,2,2,2))
plot(n, prob.alt.sg(n, null.delta1,
theta=0,
alpha=alpha.set),
type="l",
col="red",
lwd=2,
ylim=c(0,0.055),
ylab="",
main=TeX(r'(Type 1 Error $\beta_0$)'))
lines(n, prob.alt.sg(n, null.delta2,
theta=0,
alpha=alpha.set),
col="blue",
lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
theta=0,
alpha=alpha.set),
col="black",
lwd=2)
# lines(n, prob.alt.sg(n, null.delta3,
# theta=0),
# col="hotpink",
# lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
theta=0,
alpha=alpha.set),
col="hotpink",
lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
theta=0,
alpha=alpha.set),
col="orange",
lwd=2)
abline(h=0.05, col="black", lty=2)
plot(n, prob.alt.sg(n, null.delta1,
theta=1,
alpha=alpha.set),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Power $\beta_1$)'))
lines(n, prob.alt.sg(n, null.delta2,
theta=1,
alpha=alpha.set),
col="blue",
lwd=2)
lines(n, prob.alt.sg(n, null.delta0,
theta=1,
alpha=alpha.set),
col="black",
lwd=2)
# lines(n, prob.alt.sg(n, null.delta3,
# theta=1),
# col="hotpink",
# lwd=2)
lines(n, prob.alt.sg(n, null.delta4,
theta=1,
alpha=alpha.set),
col="hotpink",
lwd=2)
lines(n, prob.alt.sg(n, null.delta5,
theta=1,
alpha=alpha.set),
col="orange",
lwd=2)
legend(70, 0.8,
col=c("red",
"black",
"blue",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1),
cex=0.7,
lwd=2,
legend=c("Point Null: Constant at 0",
"Original SGPV: Constant at 0.9",
TeX(r'($n^{-1/4}$)'),
TeX(r'($n^{-1/3}$)'),
TeX(r'($n^{-1/2}$)')))
par(mfrow=c(1,2))
par(mar=c(2,2,2,2))
plot(n, prob.null.sg(n, null.delta1,
theta=0,
alpha=alpha.set),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main= TeX(r'(Probability of Null $\omega_0$)'))
lines(n, prob.null.sg(n, null.delta2,
theta=0,
alpha=alpha.set),
col="blue",
lwd=2)
lines(n, prob.null.sg(n, null.delta0,
theta=0,
alpha=alpha.set),
col="black",
lwd=2)
# lines(n, prob.null.sg(n, null.delta3,
# theta=0),
# col="hotpink",
# lwd=2)
lines(n,prob.null.sg(n, null.delta4,
theta=0,
alpha=alpha.set),
col="hotpink",
lwd=2)
lines(n, prob.null.sg(n, null.delta5,
theta=0,
alpha=alpha.set),
col="orange",
lwd=2)
plot(n, prob.null.sg(n, null.delta1,
theta=1,
alpha=alpha.set),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Probability of Null $\omega_1$)'))
lines(n, prob.null.sg(n, null.delta2,
theta=1,
alpha=alpha.set),
col="blue",
lwd=2)
lines(n, prob.null.sg(n, null.delta0,
theta=1,
alpha=alpha.set),
col="black",
lwd=2)
# lines(n, prob.null.sg(n, null.delta3,
# theta=1),
# col="hotpink",
# lwd=2)
lines(n, prob.null.sg(n, null.delta4,
theta=1,
alpha=alpha.set),
col="hotpink",
lwd=2)
lines(n, prob.null.sg(n, null.delta5,
theta=1,
alpha=alpha.set),
col="orange",
lwd=2)
legend(70, 0.75,
col=c("red",
"black",
"blue",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1),
cex=0.7,
lwd=2,
legend=c("Point Null: Constant at 0",
"Original SGPV: Constant at 0.9",
TeX(r'($n^{-1/4}$)'),
TeX(r'($n^{-1/3}$)'),
TeX(r'($n^{-1/2}$)')))
par(mfrow=c(1,2))
par(mar=c(2,2,2,2))
plot(n,prob.incon.sg(n, null.delta1,
theta=0,
alpha=alpha.set),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Probability of Inconclusive $\gamma_0$)'))
lines(n,prob.incon.sg(n, null.delta2,
theta=0,
alpha=alpha.set),
col="blue",
lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
theta=0,
alpha=alpha.set),
col="black",
lwd=2)
# lines(n,prob.incon.sg(n, null.delta3,
# theta=0),
# col="hotpink",
# lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
theta=0,
alpha=alpha.set),
col="hotpink",
lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
theta=0,
alpha=alpha.set),
col="orange",
lwd=2)
plot(n,prob.incon.sg(n, null.delta1,
theta=1,
alpha=alpha.set),
type="l",
col="red",
lwd=2,
ylim=c(0,1),
ylab="",
main=TeX(r'(Probability of Inconclusive $\gamma_1$)'))
lines(n,prob.incon.sg(n, null.delta2,
theta=1,
alpha=alpha.set),
col="blue",
lwd=2)
lines(n,prob.incon.sg(n, null.delta0,
theta=1,
alpha=alpha.set),
col="black",
lwd=2)
# lines(n,prob.incon.sg(n, null.delta3,
# theta=1),
# col="hotpink",
# lwd=2)
lines(n,prob.incon.sg(n, null.delta4,
theta=1,
alpha=alpha.set),
col="hotpink",
lwd=2)
lines(n,prob.incon.sg(n, null.delta5,
theta=1,
alpha=alpha.set),
col="orange",
lwd=2)
legend(70, 0.55,
col=c("red",
"black",
"blue",
"hotpink",
"orange"),
box.lty=0,
lty=c(1,
1,
1,
1,
1),
cex=0.7,
lwd=2,
legend=c("Point Null: Constant at 0",
"Original SGPV: Constant at 0.9",
TeX(r'($n^{-1/4}$)'),
TeX(r'($n^{-1/3}$)'),
TeX(r'($n^{-1/2}$)')))
1.4.3 COVID Clinical Trail Example
- Randomized 1,591 patients to ivermectin treatment or placebo
- Mean time spent unwell was estimated using a longitudinal ordinal regression model; range was 0 to 14 days
- Patients reported each day their symptoms and severity, health care visits, and medications.
- “The difference in the amount of time spent feeling unwell with COVID was estimated to be 0.49 days in favor of ivermectin with a 95% credible interval of (0.15, 0.82).”
#Setup
theta0 = 0
point = 0.49
data.lo = 0.15
data.hi = 0.82
n=1591
data.sd=0.17
library(sgpv)1 day difference
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -1,
null.hi = 1)## $p.delta
## [1] 1
##
## $delta.gap
## [1] NAShrinking
k=1.5
V=data.sd^2
shrink.int=(k-(k/2))*sqrt(V)/(n^(1/2))+(k/2)
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -1*shrink.int,
null.hi = shrink.int)## $p.delta
## [1] 0.9002933
##
## $delta.gap
## [1] NAk=0.75
shrink.int=(k-(k/2))*sqrt(V)/n^(1/2)+(k/2)
sgpvalue(est.lo = data.lo,
est.hi = data.hi,
null.lo = -1*shrink.int,
null.hi = shrink.int)## $p.delta
## [1] 0.3382063
##
## $delta.gap
## [1] NA